/*
 *  Pulsar Scintillation Processing Program.
 *      
 *  Copyright (C) 2012 Huang, Feng <446140689@qq.com> & Wang, Ziyou <dxwzymail@gmail.com>
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License along
 *  with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>//standard buffered input/output
#include <string.h>//string operations
#include <math.h>//mathematical declarations
#include <cpgplot.h>//header file (cpgplot.h) are provided that allow PGPLOT to be called from a C or C++ program
#include <unistd.h>//standard symbolic constants and types
#include <dirent.h>//format of directory entries
#include <sys/stat.h>//getting information about files attributes
#include <sys/types.h>//Defines primitive system data types
#include <fftw3.h>//This library (FFTW) computes Fast Fourier Transforms in one or more dimensions.
#include "ginv.c"//calculate generalized inverse of a matrix
#include "trmul.c"//matrix multiplication function
#define e 2.718281828

int nfile=-1;//the number of files that need to be handled

int main()
{
	int i,j；//Induction variable
	int Filesearch();//File search function
	char path[255];//store the path that you entered 
	void Choosefile();//Choose file function

	/*** search files ***/
	printf("This is a Pulsar Scintillation Processing Program.\nPlease enter the path of the file or folder you wish to open(.czrf format only):\n");
	while (nfile<=0)//loop until find at least one .czrf file
	{
		scanf("%s",path);//enter the path
		getchar();//Using getchar() gets the 'Enter' after input
		nfile=Filesearch(path);//call file search function
		if (nfile<=0) 
		{
			printf("find no .czrf file(s).Please input again\n");
		}
	}
	printf("find %d .czrf file(s)\n",nfile);

	/*** processing stage***/
	for (;;)	//infinite loop
	{
		printf("############Congratulations#############\n"
				"Please choose the processing procedure:\n"
				"1.Data Pre-Processing(must be run first)\n"
				"2.Dynamic Spectra\n"
				"3.Auto-Correlation Function\n"
				"4.Scintillation Parameters\n"
				"5.Secondary spectra\n"
				"6.exit\n");
		char choose[1];
		scanf("%s",choose);//choose program Functions
		getchar();//Using getchar() gets the 'Enter' after input

		if (strcmp(choose,"1") == 0) //Data Pre-Processing
		{
			FILE *fp1 = fopen("filenames.txt","r");
			if (fp1 == NULL)
			{ 
				printf("Can not open file filenames.txt,please check!\n");
				exit(1);
			}
			char choice[1];
			int algorithm=1,wrong=1;
			printf("Please choose the algorithm you want to use:\n1.40 percent of the peak\n2.sum of bins\n");
			while (wrong)
			{
				scanf("%s",choice);//choose algorithm
				getchar();//Using getchar() gets the 'Enter' after input
				if (strcmp(choice,"1")==0)
				{
					wrong=0;algorithm=1;
				}
				else if (strcmp(choice,"2")==0)
				{
					wrong=0;algorithm=2;
				}
				else
				{
					wrong=1;
					printf("Wrong input,Please choose again:");
				}
			}
			printf("Please input the Distance of the pulsar:");
			char temp[5];
			float Dist;
			scanf("%s",temp);
			Dist=atof(temp);//Convert string to double
			printf("### This process will take about %d minutes,please wait...###\n",nfile*4);//assume that each file takes four minutes
			for (i=0;i<nfile;i++)
			{
				char filename[64],filename2[64];
				char codevap[128],codepdv[128],codepam[128];
				void Convertdata(char*,int,float);//Data conversion function（filename，algorithm，distance）
				fscanf(fp1,"%s\n",filename);//Read a filename (including extension)
				//Remove File Extensions
				for (j=0;j<strlen(filename)-5;j++)
				{
					filename2[j]=filename[j];
					filename2[j+1]='\0';
				}
				//output tables of parameters and derived values
				sprintf(codevap,"vap -c \"name nsub nchan nbin bw freq length tsub dm mjd\" %s > %s.txt",filename,filename2);
				system(codevap);
				//command line general purpose data reduction
				sprintf(codepam,"pam -DpFT -e dm0pFT %s",filename);
				system(codepam);
				sprintf(codepam,"pam -Dp -t 2 -e czdm0prft2 %s",filename);
				system(codepam);
				//display Pulse widths and mean flux density;Convert data to text file format
				sprintf(codepdv,"pdv -t %s.dm0pFT >> %s.txt",filename2,filename2);
				system(codepdv);
				sprintf(codepdv,"pdv -t %s.czdm0prft2 >> %s.txt",filename2,filename2);
				system(codepdv);
				sprintf(codepdv,"pdv -f %s.czdm0prft2 > %sS.txt",filename2,filename);
				system(codepdv);
				printf("%sS.txt wrriten to disk\n",filename);	
				printf("%s.txt wrriten to disk\n",filename2);
				printf("converting data...\n");
				Convertdata(filename2,algorithm,Dist);//transforming raw data
				printf("### There remains about %d minutes ###\n",4*(nfile-i-1));
			}
			fclose(fp1);	
		}

		if (strcmp(choose,"2") == 0)//Dynamic Spectra
		{
			char filename[64];
			int DrawDynamicSpectra();//Draw Dynamic Spectra function
			Choosefile(nfile,filename);//Choosefile（number of files，filename）
			if (strcmp(filename,"ffff")==0)//Plot more than one file;batch
			{
				char input[2];
				int j,k,nx,ny,size,wrong=1;
				//Enter the number of rows
				printf("Please input the nx(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if ((atoi(input)>0 && atoi(input)<11))
					{
						wrong=0;nx=atoi(input);//atoi():Convert string to integer
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input nx again(1-10):");
					}
				}
				wrong=1;
				//Enter the number of columns
				printf("Please input the ny(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;ny=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input ny again(1-10):");
					}
				}
				wrong=1;
				//Enter the size of panel
				printf("Please input the size of panel(inch:1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;size=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input size of panel again(1-10):");
					}
				}
				FILE *fp = fopen("filenames.txt","r");
				if (fp == NULL)
				{
					printf("Can not open file %s,please check!\n","filenames.txt");
		 			return;
				}
				int n=0;
				for (i=0;i<(int)(nfile/(nx*ny))+1;i++)
				{
					if (n==nfile) break;
					cpgbeg(0,"?",nx,ny);
					cpgpap(size*nx,((float)ny)/((float)nx));
					for(j=0;j<ny;j++)
					{
						for(k=0;k<nx;k++)
						{
							if (n==nfile) break;
							fscanf(fp,"%s\n",filename);
							printf("%s\n",filename);//output raw data filename
							sprintf(filename,"%s.txt",filename);//Adding file extensions
							DrawDynamicSpectra(filename);
							cpgpanl(k+1,j+1);
							n++;
						}
					}
				}
			}
			//just one file
			else
			{
				printf("chosed file:%s\n",filename);//output raw data filename
				sprintf(filename,"%s.txt",filename);//Adding file extensions
				cpgbeg(0,"?",1,1);
				cpgpap(7.0,1.0);
				DrawDynamicSpectra(filename);
			}
			cpgend();
		}	

		if (strcmp(choose,"3") == 0) //Auto-Correlation Function
		{
			char filename[128];
			int DrawACF();//Draw ACF function
			Choosefile(nfile,filename);//Choosefile（number of files，filename）
			if (strcmp(filename,"ffff")==0)//Plot more than one file;batch
			{
				char input[2];
				int j,k,nx,ny,size,wrong=1;
				//Enter the number of rows
				printf("Please input the nx(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if ((atoi(input)>0 && atoi(input)<11))
					{
						wrong=0;nx=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input nx again(1-10):");
					}
				}
				wrong=1;
				//Enter the number of columns
				printf("Please input the ny(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;ny=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input ny again(1-10):");
					}
				}
				wrong=1;
				//Enter the size of panel
				printf("Please input the size of panel(inch:1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;size=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input size of panel again(1-10):");
					}
				}
				FILE *fp = fopen("filenames.txt","r");
				if (fp == NULL)
				{
					printf("Can not open file %s,please check!\n","filenames.txt");
		 			return;
				}
				int n=0;
				for (i=0;i<(int)(nfile/(nx*ny))+1;i++)
				{
					if (n==nfile) break;
					cpgbeg(0,"?",nx,ny);
					cpgpap(size*nx,((float)ny)/((float)nx));
					for(j=0;j<ny;j++)
					{
						for(k=0;k<nx;k++)
						{
							if (n==nfile) break;
							cpgpanl(k+1,j+1);
							fscanf(fp,"%s\n",filename);
							printf("%s\n",filename);//output raw data filename
							sprintf(filename,"%s.txt",filename);//Adding file extensions
							DrawACF(filename);
							cpgpanl(k+1,j+1);
							n++;
						}
					}
				}
			}
			//just one file
			else
			{
				printf("chosed file:%s\n",filename);//output raw data filename
				sprintf(filename,"%s.txt",filename);//Adding file extensions
				cpgbeg(0,"?",1,1);
				cpgpap(7.0,1.0);
				DrawACF(filename);
			}
			cpgend();
		}

		if (strcmp(choose,"4") == 0) //Scintillation Parameters
		{
			int DrawScinParas();//Draw Scintillation Parameters function
			DrawScinParas();
			
		}

		if (strcmp(choose,"5") == 0) //Secondary spectra
		{
			char filename[128];
			int DrawSecondaryspectra();//Draw Secondary spectra function
			Choosefile(nfile,filename);//Choose file(number of files,filename)
			if (strcmp(filename,"ffff")==0)//Plot more than one file;batch
			{
				char input[2];
				int j,k,nx,ny,size,wrong=1;
				//Enter the number of rows
				printf("Please input the nx(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if ((atoi(input)>0 && atoi(input)<11))
					{
						wrong=0;nx=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input nx again(1-10):");
					}
				}
				wrong=1;
				//Enter the number of columns
				printf("Please input the ny(1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;ny=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input ny again(1-10):");
					}
				}
				wrong=1;
				//Enter the size of panel
				printf("Please input the size of panel(inch:1-10):");
				while (wrong)//If input errors,loop to ask user input
				{
					scanf("%s",input);
					getchar();//Using getchar() gets the 'Enter' after input
					if (atoi(input)>0 && atoi(input)<11)
					{
						wrong=0;size=atoi(input);
					}
					else
					{
						wrong=1;
						printf("Wrong input,Please input size of panel again(1-10):");
					}
				}
				FILE *fp = fopen("filenames.txt","r");
				if (fp == NULL)
				{
					printf("Can not open file %s,please check!\n","filenames.txt");
		 			return;
				}
				int n=0;
				for (i=0;i<(int)(nfile/(nx*ny))+1;i++)
				{
					if (n==nfile) break;
					cpgbeg(0,"?",nx,ny);
					cpgpap(size*nx,((float)ny)/((float)nx));
					for(j=0;j<ny;j++)
					{
						for(k=0;k<nx;k++)
						{
							if (n==nfile) break;
							fscanf(fp,"%s\n",filename);
							printf("%s\n",filename);//output raw data filename
							sprintf(filename,"%s.txt",filename);//Adding file extensions(.txt)
							DrawSecondaryspectra(filename);
							cpgpanl(k+1,j+1);
							n++;
						}
					}
				}
			}
			//just one file
			else
			{
				printf("chosed file:%s\n",filename);//output raw data filename
				sprintf(filename,"%s.txt",filename);//Adding file extensions
				cpgbeg(0,"?",1,1);
				cpgpap(7.0,1.0);
				DrawSecondaryspectra(filename);
			}
			cpgend();

		}

		if (strcmp(choose,"6") == 0)	//Exit
		{
			exit(1);
		}	
	}
}


void Choosefile(int nfile,char *filename)//choose file function（number of files，filename）
{
	char choice[128];
	_Bool wrong=1;
	int i;
	FILE *fp = fopen("filenames.txt","r");
	if (fp == NULL)
	{
		 printf("Can not open file %s,please check!\n","filenames.txt");
		 return;
	}

	if (nfile>1)//More than one file
	{
		for (i=0;i<nfile;i++)
		{
			fscanf(fp,"%s\n",filename);//read the filename
			printf("%d.%s\n",i+1,filename);//show sequence number and filename
		}
		fclose(fp);
		printf("Please choose the file you want to display(0 for batch):");
		while (wrong)//Loop until user input is correct
		{
			scanf("%s",choice);//choose batch or single plot
			getchar();//Using getchar() gets the 'Enter' after input
			if ((atoi(choice)>0)&&(atoi(choice)<nfile+1))//user input is correct
			{
				wrong=0;
				FILE *fp = fopen("filenames.txt","r");
				if (fp == NULL)
				{
					printf("Can not open file %s,please check!\n","filenames.txt");
		 			return;
				}
				for (i=0;i<atoi(choice);i++)
				{
					fscanf(fp,"%s\n",filename);//read the filename
				}
				fclose(fp);
			}
			else if(strcmp(choice,"0")==0)//If user chooses 0, change filename to ffff
			{
				for (i=0;i<4;i++)
				{
					filename[i]='f';
					filename[i+1]='\0';
					wrong=0;
				}
			}
			else//If input errors,ask user input again
			{
				wrong=1;
				printf("Wrong input,Please choose again:");
			}
		}
		return;
	}
	else
	{
		fscanf(fp,"%s\n",filename);
		fclose(fp);
	}	
}

int Filesearch(char*dir)//search file (path)
{
	int n=0,i,j,k;
	char *filename,path[100];
	DIR *dp;
	struct dirent *entry;
	struct stat fileStat;
	
	if(stat(dir, &fileStat) < 0)//if not vilid path
	{ 
		printf("The path that you entered is invalid.\n");
		n=-1;
	}

	else if(S_ISDIR(fileStat.st_mode))// if a directory
	{
		if((dp = opendir(dir)) == NULL) 
		{
			fprintf(stderr,"cannot open directory: %s\n", dir);
			return(0);
		}
		chdir(dir);//changes the current directory to the specified directory
		FILE *fp = fopen("filenames.txt","w");//filenames.txt is used to record all the filenames of the specified directory
		if (fp == NULL)
		{
			 printf("Can not open file %s,please check!\n","filenames.txt");
			 return(-1);
		}
		printf("Directory scan of %s\n",dir);
		while((entry = readdir(dp)) != NULL) 
		{
			lstat(entry->d_name,&fileStat);
			if(S_ISDIR(fileStat.st_mode)) 
			{			 
				continue;
			}
			else if (strcmp(&entry->d_name[strlen(entry->d_name)-5],".czrf")==0)
			{
				printf("%d--%s\n",n+1,entry->d_name);
				fprintf(fp,"%s\n",entry->d_name);//write to the filename
				n=n+1;
			}
		}
		closedir(dp);
		fclose(fp);
	}
	else
	{
		if(strcoll(&dir[strlen(dir)-5],".czrf")==0)//if a valid filename
		{
			n=1;
			filename=strrchr(dir,'/')+1;
			//get the path of the file
			for (i=0;i<strlen(dir)-strlen(filename);i++)
			{
				path[i]=dir[i];
				path[i+1]='\0';
			}
			chdir(path);//changes the current directory to the specified directory
			FILE *fp1 = fopen("filenames.txt","w");
			if (fp1 == NULL)
			{	
				printf("Can not open file %s,please check!\n","filenames.txt");
				return(-1);
			}
			printf("%s\n",filename);
			fprintf(fp1,"%s\n",filename);
			fclose(fp1);
		}
		else
		{
			n=0;
		}
	}
	return n;//return the number of files
}

int DrawDynamicSpectra(const char *filename)//Draw Dynamic Spectra function（filename）
{
	int i,j,nsub,nchan,nbin;
	char filename11[64],name[11],title[64];//filename，pulsar name，dynamic spectrum of the title
	float bw,freq,length,tsub,dm,mjd,Dist;

	FILE *fin=fopen(filename,"r");
	if(fin==NULL)
	{
		printf("Can not open file %s,please check!\n",filename);
		exit(1);
	}
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd,&Dist);
	float a[nsub],b[nchan],A[nsub][nchan],Agray[nsub][nchan],AMAX=0;
	float tr[6]={freq-bw/2-bw/(nchan*2),bw/nchan,0,0.0-length/(nsub*120),0,length/(nsub*60)};
	for(i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
			fscanf(fin,"%f %f %f\n",&a[i],&b[j],&A[i][j]);//read data
			AMAX=(AMAX>A[i][j]?AMAX:A[i][j]);//Find the intensity of the maximum
		}
	}
	float black=0.85,white=0.01;//Set the gray scale
	for(i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
			Agray[i][j]=(A[i][j]-AMAX*white)*(50/(AMAX*black-AMAX*white))+21;//Transforming intensity values,range 21-71
		}
	}
	fclose(fin);
	cpgsch(1.6);
	cpgslw(4);
	cpgenv(freq-0.5*bw-1.5,freq+0.5*bw+1.5,-1.5,length/60.0+1.5,2,0);
	sprintf(title,"%-15s%.2f",name,mjd);//Generate dynamic spectrum of the title
	cpglab("Freq(MHz)", "Time(min)",title);
	//Set of gray-scale between 21 per cent and 71 per cent of the maximum.
	for (i=21;i<71;i++)
	{
		cpgscr(i,(i-21)/50,(i-21)/50,(i-21)/50);
	}
	cpgscir(21,70);
	cpggray(*Agray,nchan,nsub,1,nchan,1,nsub,70,21,tr);//Draw gray scale image
}

int DrawACF(const char *filename)//Draw ACF function（filename）
{
	int i,j,k,nsub,nchan,nbin;//Induction variable(i,j,k)
	char filename11[64],name[11],title[64];//filename，pulsar name，dynamic spectrum of the title
	float bw,freq,length,tsub,dm,mjd,Dist,MIN(float,float),MAX(float,float),ACF();

	FILE *fin=fopen(filename,"r");
	if(fin==NULL)
	{
		printf("Can not open file %s,please check!\n",filename);
		exit(1);
	}
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd,&Dist);
	float a[nsub],b[nchan],A[nsub][nchan],Agray[nsub][nchan],S[nsub][nchan],AMAX=0,AVERAGE=0,td0,vd0;
	int zeron=0;//the number of intensity value of 0
	for(i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
			fscanf(fin,"%f %f %f\n",&a[i],&b[j],&A[i][j]);
			AMAX=(AMAX>A[i][j]?AMAX:A[i][j]);//Find the intensity of the maximum
			AVERAGE+=A[i][j];
			if (A[i][j]==0)
			{
				zeron++;
			}
		}
	}
	AVERAGE=AVERAGE/(nsub*nchan-zeron);//The average intensity values
	//normalize intensity values 
	for (i=0;i<nsub;i++)
	{
		for (j=0;j<nchan;j++)
		{
			if(A[i][j]==0)
			{
				S[i][j]=0.0;
			}
			else S[i][j]=A[i][j]-AVERAGE;
		}
	}
	fclose(fin);
	float C[nsub][nchan*2-1],xr[nchan*2-1],kx[nchan*2-1],yr[nsub],ky[nsub],xrMIN=1000,yrMIN=1000;
	//normalized two-dimensional ACFs values,xr is x-axis range in the bottom of figure,kx is y-axis range in the bottom of figure,yr is x-axis range in the right of figure,ky is y-axis range in the right of figure
	for (i=0;i<nsub;i++)
	{
		for (j=-nchan+1;j<nchan;j++)
		{
			C[i][j+nchan-1]=ACF(j,i,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);//normalize ACF values 
		}
	}
	for (j=0;j<nchan;j++)
	{
		kx[j]=(j-nchan+1)/(nchan/bw);
		kx[j+nchan-1]=j/(nchan/bw);
		xr[j+nchan-1]=ACF(j,0,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		xr[j]=ACF(nchan-j-1,0,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		xrMIN=MIN(xr[j],xrMIN);
	}
	for (j=0;j<nchan;j++)//Find half-width at half-maximum of the ACF along the frequency-lag axis at zero time lag
	{
		if (xr[j+nchan]<0.5)
		{
			vd0=j*(bw/nchan);
			break;
		}
	}
	for (i=0;i<nsub;i++)
	{
		ky[i]=i*(length/(nsub*60));
		yr[i]=ACF(0,i,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		yrMIN=MIN(yr[i],yrMIN);
	}
	for (i=0;i<nsub;i++)//Find the time lag at zero frequency lag at the point where the ACF is 1/e of the maximum
	{
		if (yr[i]<1/e)
		{
			td0=i*(length/(nsub*60));
			break;
		}
	}

	printf("td0:%f\tvd0:%f\n",td0,vd0);

	int t=td0/2,v=vd0;//Set the fitting range
	if (t<4) t=4;//In order to avoid fitting too few data
	if (v<1) v=2;
	double Clog[t*2*v],D1[t*2*v][3],DD[3][t*2*v],U[t*2*v][t*2*v],V[3][3],X[3],Clog2[t*2*v];
	for (i=0,k=0;i<t;i++)
	{
		for (j=-v;j<v;j++)
		{
			if (C[i][j+nchan]<=0) Clog[k]=log(0.000001);
			else Clog[k]=log(C[i][j+nchan]);
			D1[k][0]=(bw/nchan*j)*(bw/nchan*j);
			D1[k][1]=(bw/nchan*j)*(i*(length/(nsub*60)));
			D1[k][2]=(i*length/(nsub*60))*(i*length/(nsub*60));
			k++;
		}
	}
	
	ginv(D1,t*2*v,3,DD,0.0001,U,V,t*2*v+1);//call function of calculate generalized inverse of a matrix 
	trmul(DD,Clog,3,t*2*v,1,X);//call matrix multiplication function
	trmul(D1,X,t*2*v,3,1,Clog2);

	double SD=0.0;//fitting error
	for (i=0;i<t*2*v;i++)
	{
		SD+=(Clog2[i]-Clog[i])*(Clog2[i]-Clog[i]);
	}
	SD=sqrt(SD/(t*2*v));
	printf("c1:%f\tc2:%f\tc3:%f\tSD:%f\n",-X[0],-X[1],-X[2],SD);
	printf("td:%f\tSD:%f\tvd:%f\tSD:%f\n",sqrt(-1/(X[2])),-0.5*pow(-X[2],-1.5)*SD,sqrt(-log(2)/(X[0])),-0.5*pow(-X[0],-1.5)/(sqrt(log(2)))*SD);

	float gsxr[nchan*2-1],gsyr[nsub];
	for (j=0;j<nchan;j++)
	{
		gsxr[j+nchan-1]=exp((X[0])*(bw/nchan*j)*(bw/nchan*j));
		gsxr[j]=exp((X[0])*(bw/nchan*(j-nchan+1))*(bw/nchan*(j-nchan+1)));
	}
	for (i=0;i<nsub;i++)
	{
		gsyr[i]=exp((X[2])*(i*length/(nsub*60))*(i*length/(nsub*60)));
	}

	//draw two-dimensional ACFs of the dynamic spectra
	float tr[6]={-bw-bw/(nchan*2),bw/nchan,0,0.0-length/(nsub*120),0,length/(nsub*60)};
	float ALEV[21]={-1,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1};
	cpgslw(3);
	cpgsls(1);
	cpgsch(1);
	cpgsvp(0.1,0.7,0.3,0.9);
	cpgswin(-bw-1.5,bw+1.5,-1.5,length/60.0+1.5);
	cpgbox("bcst",0.,10,"bcnst",0.0,0.1);
	sprintf(title,"%-15s %.2f",name,mjd);
	cpglab("","Time lag(min)",title);
	cpgcont(*C,nchan*2-1,nsub,1,nchan*2-1,1,nsub,ALEV,21,tr);
	 
	//one-dimensional ACFs at zero time lag,at the bottom
	cpgslw(3);
	cpgsch(1);
	cpgsvp(0.1,0.7,0.1,0.25);
	cpgswin(-bw-1.5,bw+1.5,xrMIN-0.3,1.3);
	cpgbox("bcnst",0.,10,"bcnst",0.,0.);
	cpglab("Freq lag(MHz)","","");
	cpgline(nchan*2-1,kx,xr);
	cpgsls(4);
	float kx2[nchan-1],gsxr2[nchan-1];
	for (i=0;i<nchan-1;i++)
	{
		kx2[i]=kx[2*i];
		gsxr2[i]=gsxr[2*i];
	}
	cpgslw(4);
	cpgpt(nchan-1,kx2,gsxr2,-1);
		
	//one-dimensional ACFs at zero frequency lag,at the right
	cpgslw(4);
	cpgsch(1);
	cpgsls(1);
	cpgsvp(0.75,0.9,0.3,0.9);
	cpgswin(yrMIN-0.05,1.,-1.5,length/60.0+1.5);
	cpgbox("bcnst",0.,10.,"bcst",0.,0.);
	cpgline(nsub,yr,ky);
	cpgsls(4);
	cpgslw(4);
	cpgpt(nsub,gsyr,ky,-1);
}


int DrawSecondaryspectra(char *filename)//Draw Secondary spectra function（filename）
{
	int i,j,nsub,nchan,nbin;
	char filename11[64],name[11],title[64];
	float bw,freq,length,tsub,dm,mjd,Dist,MAX(float,float);

	FILE *fin=fopen(filename,"r");
	if(fin==NULL)
	{
		printf("Can not open file %s,please check!\n",filename);
		exit(1);
	}
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd,&Dist);
	float a[nsub],b[nchan],A[nsub][nchan],Agray[nsub][nchan],AMAX=0;
	float tr[6]={freq-bw/2-bw/(nchan*2),bw/nchan,0,0.0-length/(nsub*120),0,length/(nsub*60)};
	for(i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
		fscanf(fin,"%f %f %f\n",&a[i],&b[j],&A[i][j]);
		AMAX=(AMAX>A[i][j]?AMAX:A[i][j]);	
		}
	}

	float D[nsub][nchan],Dgray[nsub][nchan],DMAX=0.;
	// Fast Fourier Transforms in two dimensions.
	const ptrdiff_t N0 = nsub, N1 = nchan;
	fftw_plan plan;
	fftw_complex *data;
	data = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * N0 * N1);
 
	// create plan for forward DFT 
	plan = fftw_plan_dft_2d(N0, N1, data, data, FFTW_FORWARD, FFTW_ESTIMATE);
 
	// initialize data to some function my_function(x,y) 
  	double pdata=0;
	for (i = 0; i < N0; ++i)
		{
		for (j = 0; j < N1; ++j)
			{
				data[i*N1 + j][0]=A[i][j]; 
				data[i*N1 + j][1]=0;
				pdata+=data[i*N1 + j][0]*data[i*N1 + j][0]+data[i*N1 + j][1]*data[i*N1 + j][1];
	  		}
		}
	//printf("power of original data is %f\n", pdata);
 
	// compute transforms, in-place, as many times as desired 
	fftw_execute(plan);
 
	double normalization=sqrt((double)N0*N1);
	double ptransform = 0;
 
	//normalize data and calculate power of transform 
	for (i = 0; i < N0; ++i)
	{
		for (j = 0; j < N1; ++j)
		{
			data[i*N1+j][0]/=normalization;
			data[i*N1+j][1]/=normalization;
			ptransform+=data[i*N1 + j][0]*data[i*N1 + j][0]+data[i*N1 + j][1]*data[i*N1 + j][1];
			D[i][j] = sqrt(data[i*N1 + j][0]*data[i*N1 + j][0]+data[i*N1 + j][1]*data[i*N1 + j][1]);
			DMAX=MAX(D[i][j],DMAX);
		}
	}
	//printf("power of transform is %f\n", ptransform);
	
	for (i = 0; i < N0; ++i)
	{
		for (j = 0; j < N1; ++j)
		{
			if(i<N0/2) Dgray[i][j]=D[i+N0/2][j];
			else Dgray[i][j]=D[i-N0/2][j];
		}
	}
	for (i = 0; i < N0; ++i)
	{
		for (j = 0; j < N1; ++j)
		{
			if(j<N1/2) Dgray[i][j]=Dgray[i][j+N1/2];
			else Dgray[i][j]=Dgray[i][j-N1/2];
		}
	}
				
	fftw_destroy_plan(plan);
	fftw_free(data); 

	sprintf(title,"%s %f",name,mjd);
	// Use PGENV to specify the range of the axes and to draw a box, and PGLAB to label it.
	cpgenv(-0.01,0.21,-0.12,0.12,0,0);
	cpglab("fv(cycles/MHz)","ft(cycles/min)",title);
	float TR[6]={0,0.2/nchan,0,-0.11,0,0.22/nsub};
	cpggray(*Dgray,nchan,nsub,1,nchan,1,nsub,DMAX,0,TR);
}

void Convertdata(char *filename,int algorithm,float Dist)//Data conversion function（filename，algorithm，distance）
{
	int i,j,k,l,r,m,n,nsub,nchan,nbin;
	char filename1[9],filename11[40],filename2[40],name1[5],name[11],nsub1[5],nchan1[6],nbin1[5],
			bw1[3],freq1[6],length1[7],tsub1[5],dm1[3],mjd1[4],rubbish1[8],Dist1[5]="Dist";
	float bw,freq,length,tsub,dm,mjd;
	float MAX(float,float);
	sprintf(filename2,"%s.czrf.txt",filename);//Adding file extensions(.czrf.txt)
	sprintf(filename,"%s.txt",filename);//Adding file extensions(.txt)
	FILE *fin = fopen(filename,"r");
	if (fin == NULL)
	{
		 printf("Can not open file %s,please check!\n",filename);
		 exit(1);
	}
	 
	fscanf(fin,"%s %s %s %s %s %s %s %s %s %s %s\n",filename1,name1,nsub1,nchan1,nbin1,bw1,freq1,length1,tsub1,dm1,mjd1);
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd);
	fscanf(fin,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s",rubbish1,rubbish1,rubbish1,rubbish1,
			rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1);
	nsub=nsub/2;
	float a[nbin],b[nbin],c[nbin],A[nsub][nchan],bMAX=0;

	//assign initial values
	for (i=0;i<nsub;i++)
		{
			for (j=0;j<nchan;j++)
			{			 
				A[i][j]=0;
			}
		}

	for (i=0;i<nbin;i++)
		{
			fscanf(fin,"%f %f %f %f",&a[i],&a[i],&a[i],&b[i]);
			bMAX=MAX(b[i],bMAX);
		}
	fscanf(fin,"%s %s %s %s %s %s %s %s %s %s %s %s %s %s",rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,
			rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1,rubbish1);
	FILE *fpp=fopen(filename2,"w");
	if (fpp == NULL)
		{
			printf("Can not open file %s,please check!\n",filename2);
			exit(1);
		}

	fprintf(fpp,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,nsub,nchan,nbin,bw,freq,length,tsub,dm,mjd,Dist);

	if (algorithm==1)//Pulse intensities were determined by summing the data within phase ±0.4 of the predicted phase of the pulse peak 
	{
		for (i=0;i<nbin;i++)
		{
			if (b[i]==bMAX)
			{
				for (j=0;j<nbin;j++)
				{
					k=i-j;
					if(k<0) k+=nbin;
					if(b[k]<bMAX*0.4)
					{	l=k;
						if(l>i) l=nbin-l+i;
						else l=i-l;
						break;
					}
				}
				for (j=0;j<nbin;j++)
				{	k=i+j;
					if(k>nbin-1) k-=nbin;
					if(b[k]<bMAX*0.4)
					{	r=k;
						if(r<i) r=nbin-i+r;
						else r=r-i;
						break;
					}
				}
				break;
			}
			}

		for(k=0;k<nsub;k++)
		{
		for (j=0;j<nchan;j++)
		{	bMAX=0.;
			for (i=0;i<nbin;i++)
			{
				fscanf(fin,"%f %f %f %f",&a[i],&a[i],&a[i],&b[i]);
				bMAX=MAX(b[i],bMAX);
			}
			for (i=0;i<nbin;i++)
			{
				if (b[i]==bMAX) break; 
			}
			for(m=0;m<l;m++)
			{	n=i-m;
				if(n<0) n+=nbin;
				A[k][j]=A[k][j]+b[n];
			}
			for(m=0;m<r;m++)
			{	n=i+m;
				if(n>nbin-1) n-=nbin;
				A[k][j]=A[k][j]+b[n];
			}
			A[k][j]=A[k][j]-bMAX;
			if(A[k][j]<0) A[k][j]=0.001;
			fprintf(fpp,"%d %d %f\n",k,j,A[k][j]);
		}
		}
	}

	if (algorithm==2)//summing all the values
	{
		for(k=0;k<nsub;k++)
		{
			for (j=0;j<nchan;j++)
			{
				for (i=0;i<nbin;i++)
				{
					fscanf(fin,"%f %f %f %f",&a[i],&a[i],&a[i],&b[i]);
					A[k][j]+=b[i];
				}
				if(A[k][j]<0) A[k][j]=0.001;
				fprintf(fpp,"%d %d %f\n",k,j,A[k][j]);
			}
		}
	}

	if (algorithm==3)
	{
		/*for(k=0;k<nsub;k++)
		{
			for (j=0;j<nchan;j++)
			{
				for (i=0;i<nbin;i++)
				{
					fscanf(fin,"%f %f %f %f",&a[i],&a[i],&a[i],&b[i]);
					binmax=MAX(binmax,b[i]);
				}
				for (i=0;i<nbin;i++)
				{
					if (binmax==b[i])
					{maxn=i;break;}
				}
				for (i=maxn;i<nbin;i++)
				{
					if (binmax==b[i];
					maxn=i;
				}
				if(A[k][j]<0) A[k][j]=0.001;
				fprintf(fpp,"%d %d %f\n",k,j,A[k][j]);
			}
		}*/
	}

	fclose(fin);
	fclose(fpp);
}

int DrawScinParas()//Draw Scintillation Parameters function
{
	int i,j,Gsacf(char*,float*,float*,float*,float*,float*,char*);
	char filename[64],filename2[64],title[32],name[11];
	float 	tdy[nfile],td[nfile][2],tdmin=100,tdmax=0.,tdy1[nfile+1],tdy2[nfile+1],tdsum=0.0,tdave,tdSDs,tdSDave=0.0,
			vdy[nfile],vd[nfile][2],vdmin=100,vdmax=0.,vdy1[nfile+1],vdy2[nfile+1],vdsum=0.0,vdave,vdSDs,vdSDave=0.0,
			dtdvy[nfile],dtdv[nfile][2],dtdvmin=100,dtdvmax=0.,dtdvy1[nfile+1],dtdvy2[nfile+1],dtdvsum=0.0,dtdvave,dtdvSDs,dtdvSDave=0.0,
			ry[nfile],r[nfile][2],rmin=100,rmax=0.,ry1[nfile+1],ry2[nfile+1],rsum=0.0,rave,rSDs,rSDave=0.0,
			uy[nfile],u[nfile][2],umin=100,umax=0.,uy1[nfile+1],uy2[nfile+1],usum=0.0,uave,uSDs,uSDave=0.0,
			try[nfile],tr[nfile][2],trmin=100,trmax=0.,try1[nfile+1],try2[nfile+1],trsum=0.0,trave,trSDs,trSDave=0.0,
			logcn2y[nfile],logcn2[nfile][2],logcn2min=100,logcn2max=0.,logcn2y1[nfile+1],logcn2y2[nfile+1],logcn2sum=0.0,logcn2ave,logcn2SDs,logcn2SDave=0.0,
			Vsy[nfile],Vs[nfile][2],Vsmin=100,Vsmax=0.,Vsy1[nfile+1],Vsy2[nfile+1],Vssum=0.0,Vsave,VsSDs,VsSDave=0.0,
			Sy[nfile],Smin=10000,Smax=0.,Ssum=0.0,Save=0.0,
			v,Dist,MJD[nfile],c[3],
			MIN(float,float),MAX(float,float),xr[nfile+1],y,SD[nfile],ReadS(char*,char*);
	FILE *fp1 = fopen("filenames.txt","r");
	if (fp1 == NULL)
	{ 
		printf("Can not open file filenames.txt,please check!\n");
		exit(1);
	}
	printf("NO.\tMJD\t\ttd\tvd\tdtdv\t|r|\tu\ttr\tlogcn2\tVs\tS\n");
	for (i=0;i<nfile;i++)
	{
		fscanf(fp1,"%s\n",filename);
		sprintf(filename2,"%sS.txt",filename);//Adding file extensions(S.txt)
		sprintf(filename,"%s.txt",filename);//Adding file extensions(.txt)
		Gsacf(filename,c,&SD[i],&v,&Dist,&MJD[i],name);
		Sy[i]=ReadS(filename,filename2);//Read flux den­sity
		Ssum+=Sy[i];
		Smin=MIN(Sy[i],Smin);
		Smax=MAX(Sy[i],Smax);

		//calculate Scintillation Parameters and fitting errors
		td[i][0]=sqrt(1/(c[2]));
		tdy[i]=td[i][0];
		td[i][1]=-0.5*pow(c[2],-1.5)*SD[i];
		tdmin=MIN(td[i][0],tdmin);
		tdmax=MAX(td[i][0],tdmax);
		tdsum+=td[i][0];
		tdSDave+=td[i][1];
		tdy1[i]=td[i][0]-0.5*td[i][1];
		tdy2[i]=td[i][0]+0.5*td[i][1];

		vd[i][0]=sqrt(fabs(-log(2)/(c[0])));
		vdy[i]=vd[i][0];
		vd[i][1]=-0.5*pow(c[0],-1.5)/(sqrt(log(2)))*SD[i];
		vdmin=MIN(vd[i][0],vdmin);
		vdmax=MAX(vd[i][0],vdmax);
		vdSDave+=vd[i][1];
		vdsum+=vd[i][0];
		vdy1[i]=vd[i][0]-0.5*vd[i][1];
		vdy2[i]=vd[i][0]+0.5*vd[i][1];

		dtdv[i][0]=-(c[1]/(2*(c[2])));
		dtdvy[i]=dtdv[i][0];
		dtdv[i][1]=sqrt(SD[i]*SD[i]/(4*c[2]*c[2])+c[1]*c[1]/4*pow(-c[2],-4)*SD[i]*SD[i]);
		dtdvmin=MIN(dtdv[i][0],dtdvmin);
		dtdvmax=MAX(dtdv[i][0],dtdvmax);
		dtdvsum+=dtdv[i][0];
		dtdvSDave+=dtdv[i][1];
		dtdvy1[i]=dtdv[i][0]-0.5*dtdv[i][1];
		dtdvy2[i]=dtdv[i][0]+0.5*dtdv[i][1];

		r[i][0]=fabs(c[1]/(sqrt(4*fabs((c[0])*(c[2])))));
		ry[i]=r[i][0];
		r[i][1]=sqrt((pow(c[1]*c[2]/(4*pow(c[0]*c[3],1.5)),2)+pow(0.5*pow(c[0]*c[2],-0.5),2)+pow(-0.25*c[0]*c[1]/pow(c[0]*c[2],1.5),2))*SD[i]*SD[i]);
		rmin=MIN(r[i][0],rmin);
		rmax=MAX(r[i][0],rmax);
		rsum+=r[i][0];
		rSDave+=r[i][1];
		ry1[i]=r[i][0]-0.5*r[i][1];
		ry2[i]=r[i][0]+0.5*r[i][1];

		u[i][0]=sqrt(2*v/vd[i][0]);
		uy[i]=u[i][0];
		u[i][1]=sqrt(2*v/log(2))*0.5*pow(c[0],-0.5)*SD[i];
		umin=MIN(u[i][0],umin);
		umax=MAX(u[i][0],umax);
		usum+=u[i][0];
		uSDave+=u[i][1];
		uy1[i]=u[i][0]-0.5*u[i][1];
		uy2[i]=u[i][0]+0.5*u[i][1];

		tr[i][0]=2*v*td[i][0]/vd[i][0]/24/60;
		try[i]=tr[i][0];
		tr[i][1]=2*v*SD[i]/pow(log(2),2)*sqrt(1/c[2]*c[0]*c[0]+pow(c[0],4)*pow(c[2],-3));
		trmin=MIN(tr[i][0],trmin);
		trmax=MAX(tr[i][0],trmax);
		trsum+=tr[i][0];
		trSDave+=tr[i][1];
		try1[i]=tr[i][0]-0.5*tr[i][1];
		try2[i]=tr[i][0]+0.5*tr[i][1];

		logcn2[i][0]=log(0.002*pow(v/1000,11.0/3)*pow(Dist,-11.0/6)*pow(vd[i][0],-5.0/6));
		logcn2y[i]=logcn2[i][0];
		logcn2[i][1]=1/(0.002*pow(v/1000,11.0/3)*pow(Dist,-11.0/6)*pow(vd[i][0],-5.0/6))  *0.002*pow(v/1000,11.0/3)*pow(Dist,-11.0/6)*pow(log(2),-5/12)*5/12*pow(c[0],-7/12);
		logcn2min=MIN(logcn2[i][0],logcn2min);
		logcn2max=MAX(logcn2[i][0],logcn2max);
		logcn2sum+=logcn2[i][0];
		logcn2SDave+=logcn2[i][1];
		logcn2y1[i]=logcn2[i][0]-0.5*logcn2[i][1];
		logcn2y2[i]=logcn2[i][0]+0.5*logcn2[i][1];

		Vs[i][0]=3.85*10000*sqrt((vd[i][0])*Dist)/(v/1000*(td[i][0])*60);
		Vsy[i]=Vs[i][0];
		Vs[i][1]=3.85*10000*sqrt(Dist)/v/1000 *sqrt(0.5*0.5/td[i][0]*60/vd[i][0]*vd[i][1]*vd[i][1]+0.5*0.5*vd[i][0]*pow(td[i][0]*60,-3)*td[i][1]*60*td[i][1]*60);
		Vsmin=MIN(Vs[i][0],Vsmin);
		Vsmax=MAX(Vs[i][0],Vsmax);
		Vssum+=Vs[i][0];
		VsSDave+=Vs[i][1];
		Vsy1[i]=Vs[i][0]-0.5*Vs[i][1];
		Vsy2[i]=Vs[i][0]+0.5*Vs[i][1];
		printf("%d\t%.2f\t%.2f\t%.2f\t%.3f\t%.4f\t%.1f\t%.2f\t%.2f\t%.1f\t%.1f\n",i+1,MJD[i],tdy[i],vdy[i],dtdvy[i],ry[i],uy[i],try[i],logcn2y[i],Vsy[i],Sy[i]);
	}
	
	//Sorting Scintillation Parameters in accordance with MJD from small to large
	for(i=0;i<nfile-1;i++)
	{
		for(j=nfile-1;j>=i;j--)
			{
				if (MJD[j]<MJD[j-1])
				{
					y=MJD[j];MJD[j]=MJD[j-1];MJD[j-1]=y;

					y=td[j][0];td[j][0]=td[j-1][0];td[j-1][0]=y;
					y=td[j][1];td[j][1]=td[j-1][1];td[j-1][1]=y;
					y=tdy[j];tdy[j]=tdy[j-1];tdy[j-1]=y;
					y=tdy1[j];tdy1[j]=tdy1[j-1];tdy1[j-1]=y;
					y=tdy2[j];tdy2[j]=tdy2[j-1];tdy2[j-1]=y;

					y=vd[j][0];vd[j][0]=vd[j-1][0];vd[j-1][0]=y;
					y=vd[j][1];vd[j][1]=vd[j-1][1];vd[j-1][1]=y;
					y=vdy[j];vdy[j]=vdy[j-1];vdy[j-1]=y;
					y=vdy1[j];vdy1[j]=vdy1[j-1];vdy1[j-1]=y;
					y=vdy2[j];vdy2[j]=vdy2[j-1];vdy2[j-1]=y;

					y=dtdv[j][0];dtdv[j][0]=dtdv[j-1][0];dtdv[j-1][0]=y;
					y=dtdv[j][1];dtdv[j][1]=dtdv[j-1][1];dtdv[j-1][1]=y;
					y=dtdvy[j];dtdvy[j]=dtdvy[j-1];dtdvy[j-1]=y;
					y=dtdvy1[j];dtdvy1[j]=dtdvy1[j-1];dtdvy1[j-1]=y;
					y=dtdvy2[j];dtdvy2[j]=dtdvy2[j-1];dtdvy2[j-1]=y;
				
					y=r[j][0];r[j][0]=r[j-1][0];r[j-1][0]=y;
					y=r[j][1];r[j][1]=r[j-1][1];r[j-1][1]=y;
					y=ry[j];ry[j]=ry[j-1];ry[j-1]=y;
					y=ry1[j];ry1[j]=ry1[j-1];ry1[j-1]=y;
					y=ry2[j];ry2[j]=ry2[j-1];ry2[j-1]=y;

					y=u[j][0];u[j][0]=u[j-1][0];u[j-1][0]=y;
					y=u[j][1];u[j][1]=u[j-1][1];u[j-1][1]=y;
					y=uy[j];uy[j]=uy[j-1];uy[j-1]=y;
					y=uy1[j];uy1[j]=uy1[j-1];uy1[j-1]=y;
					y=uy2[j];uy2[j]=uy2[j-1];uy2[j-1]=y;

					y=tr[j][0];tr[j][0]=tr[j-1][0];tr[j-1][0]=y;
					y=tr[j][1];tr[j][1]=tr[j-1][1];tr[j-1][1]=y;
					y=try[j];try[j]=try[j-1];try[j-1]=y;
					y=try1[j];try1[j]=try1[j-1];try1[j-1]=y;
					y=try2[j];try2[j]=try2[j-1];try2[j-1]=y;

					y=logcn2[j][0];logcn2[j][0]=logcn2[j-1][0];logcn2[j-1][0]=y;
					y=logcn2[j][1];logcn2[j][1]=logcn2[j-1][1];logcn2[j-1][1]=y;
					y=logcn2y[j];logcn2y[j]=logcn2y[j-1];logcn2y[j-1]=y;
					y=logcn2y1[j];logcn2y1[j]=logcn2y1[j-1];logcn2y1[j-1]=y;
					y=logcn2y2[j];logcn2y2[j]=logcn2y2[j-1];logcn2y2[j-1]=y;

					y=Vs[j][0];Vs[j][0]=Vs[j-1][0];Vs[j-1][0]=y;
					y=Vs[j][1];Vs[j][1]=Vs[j-1][1];Vs[j-1][1]=y;
					y=Vsy[j];Vsy[j]=Vsy[j-1];Vsy[j-1]=y;
					y=Vsy1[j];Vsy1[j]=Vsy1[j-1];Vsy1[j-1]=y;
					y=Vsy2[j];Vsy2[j]=Vsy2[j-1];Vsy2[j-1]=y;

					y=Sy[j];Sy[j]=Sy[j-1];Sy[j-1]=y;
				}
			}
	}
	
	//calculate RMS(root mean square) error
	tdave=tdsum/nfile;
	vdave=vdsum/nfile;
	dtdvave=dtdvsum/nfile;
	rave=rsum/nfile;
	uave=usum/nfile;
	trave=trsum/nfile;
	logcn2ave=logcn2sum/nfile;
	Vsave=Vssum/nfile;
	Save=Ssum/nfile;

	tdSDave=tdSDave/nfile;
	vdSDave=vdSDave/nfile;
	dtdvSDave=dtdvSDave/nfile;
	rSDave=rSDave/nfile;
	uSDave=uSDave/nfile;
	trSDave=trSDave/nfile;
	logcn2SDave=logcn2SDave/nfile;
	VsSDave=VsSDave/nfile;

	for(i=0;i<nfile;i++)
	{
		xr[i]=MJD[i]-MJD[0];
		tdSDs+=(td[i][0]-tdave)*(td[i][0]-tdave);
		vdSDs+=(vd[i][0]-vdave)*(vd[i][0]-vdave);
		dtdvSDs+=(dtdv[i][0]-dtdvave)*(dtdv[i][0]-dtdvave);
		rSDs+=(r[i][0]-rave)*(r[i][0]-rave);
		uSDs+=(u[i][0]-uave)*(u[i][0]-uave);
		trSDs+=(tr[i][0]-trave)*(tr[i][0]-trave);
		logcn2SDs+=(logcn2[i][0]-logcn2ave)*(logcn2[i][0]-logcn2ave);
		VsSDs+=(Vs[i][0]-Vsave)*(Vs[i][0]-Vsave);
	}
	tdSDs=sqrt(tdSDs/nfile);
	vdSDs=sqrt(vdSDs/nfile);
	dtdvSDs=sqrt(dtdvSDs/nfile);
	rSDs=sqrt(rSDs/nfile);
	uSDs=sqrt(uSDs/nfile);
	trSDs=sqrt(trSDs/nfile);
	logcn2SDs=sqrt(logcn2SDs/nfile);
	VsSDs=sqrt(VsSDs/nfile);

	fclose(fp1);
	
	//output Scintillation Parameters to parameters.txt
	FILE *fp2 = fopen("parameters.txt","w");
	if (fp2 == NULL)
	{ 
		printf("Can not open file filenames.txt,please check!\n");
		exit(1);
	}
	fprintf(fp2,"MJD\t\t\t\ttd\t\t\tSD\t\t\tvd\t\t\tSD\t\t\tdtdv\t\tSD\t\t\t|r|\t\t\tSD\t\t\tu\t\t\tSD\t\t\ttr\t\t\tSD\t\t\tlogcn2\t\tSD\t\t\tVs\t\t\tSD\t\t\tS\n");
	for (i=0;i<nfile;i++)
	{
		fprintf(fp2,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
				MJD[i],td[i][0],td[i][1],vd[i][0],vd[i][1],dtdv[i][0],dtdv[i][1],r[i][0],
				r[i][1],u[i][0],u[i][1],tr[i][0],tr[i][1],logcn2[i][0],logcn2[i][1],Vs[i][0],Vs[i][1],Sy[i]);
	}
	fprintf(fp2,"AVERAGED:\n");
	fprintf(fp2,"name\t\ttd\t\t\tSD\t\t\tvd\t\t\tSD\t\t\tdtdv\t\tSD\t\t\t|r|\t\t\tSD\t\t\tu\t\t\tSD\t\t\ttr\t\t\tSD\t\t\tlogcn2\t\tSD\t\t\tVs\t\t\tSD\t\t\tS\n");
	fprintf(fp2,"%s\t\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",
				name,tdave,tdSDave,vdave,vdSDave,dtdvave,dtdvSDave,rave,
				rSDave,uave,uSDave,trave,trSDave,logcn2ave,logcn2SDave,Vsave,VsSDave,Save);
	fclose(fp2);
	printf("AVERAGED:\n");
	printf("%s\t\t%.2f\t%.2f\t%.3f\t%.3f\t%.1f\t%.2f\t%.2f\t%.1f\t%.1f\n",name,tdave,vdave,dtdvave,rave,uave,trave,logcn2ave,Vsave,Save);
	xr[nfile]=-xr[nfile-1]*0.025;//abscissa of error bar in the mean value

	cpgbeg(0,"?",1,1);
	cpgpap(10.0,6.0/10.0);
	cpgslw(4);
	cpgsch(1);
	cpgsvp(0.1,0.95,0.65,0.9);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,tdmin-(tdmax-tdmin)*0.2,tdmax+(tdmax-tdmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	sprintf(title,"%s",name);
	cpglab("","td(min)",title);
	cpgslw(7);
	cpgpt(nfile,xr,tdy,-1);
	cpgslw(2);
	cpgline(nfile,xr,tdy);
	tdy1[nfile]=tdmax-tdSDs;
	tdy2[nfile]=tdmax;
	//cpgerry(nfile+1,xr,tdy1,tdy2,1.0);
	float xre[1],tdy1e[1],tdy2e[1],vdy1e[1],vdy2e[1];
	xre[0]=xr[nfile];tdy1e[0]=tdy1[nfile];tdy2e[0]=tdy2[nfile];
	cpgerry(1,xre,tdy1e,tdy2e,1.0);//draw the error bar in the mean value

	cpgslw(4);
	cpgsvp(0.1,0.95,0.35,0.60);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,vdmin-(vdmax-vdmin)*0.2,vdmax+(vdmax-vdmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","vd(MHz)","");
	cpgslw(7);
	cpgpt(nfile,xr,vdy,-1);
	cpgslw(2);
	cpgline(nfile,xr,vdy);
	vdy1[nfile]=vdmax-vdSDs;
	vdy2[nfile]=vdmax;
	vdy1e[0]=vdy1[nfile];vdy2e[0]=vdy2[nfile];
	//cpgerry(nfile+1,xr,vdy1,vdy2,1.0);
	cpgerry(1,xre,vdy1e,vdy2e,1.0);//draw the error bar in the mean value

	cpgslw(4);
	cpgsvp(0.1,0.95,0.05,0.30);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,dtdvmin-(dtdvmax-dtdvmin)*0.2,dtdvmax+(dtdvmax-dtdvmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","Drift Rate(min/MHz)","");
	cpgslw(7);
	cpgpt(nfile,xr,dtdvy,-1);
	cpgslw(2);
	cpgline(nfile,xr,dtdvy);
	float xr0[2],dtdvy0[2];
	xr0[0]=-xr[nfile-1]*0.05,xr0[1]=xr[nfile-1]+xr[nfile-1]*0.05,dtdvy0[0]=0.0,dtdvy0[1]=0.0;
	cpgsls(4);
	cpgline(2,xr0,dtdvy0);
	cpgsls(1);
	//cpgerry(nfile,xr,dtdvy1,dtdvy2,1.0);
	cpgend();

	cpgbeg(0,"?",1,1);
	cpgpap(10.0,6.0/10.0);
	cpgslw(4);
	cpgsvp(0.1,0.95,0.7,0.95);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,rmin-(rmax-rmin)*0.2,rmax+(rmax-rmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","Slope Visilbility","");
	cpgslw(7);
	cpgpt(nfile,xr,ry,-1);
	cpgslw(2);
	cpgline(nfile,xr,ry);
	//cpgerry(nfile,xr,ry1,ry2,1.0);
	
	cpgslw(4);
	cpgsvp(0.1,0.95,0.4,0.65);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,umin-(umax-umin)*0.2,umax+(umax-umin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","u","");
	cpgslw(7);
	cpgpt(nfile,xr,uy,-1);
	cpgslw(2);
	cpgline(nfile,xr,uy);
	//cpgerry(nfile,xr,uy1,uy2,1.0);

	cpgslw(4);
	cpgsvp(0.1,0.95,0.1,0.35);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,trmin-(trmax-trmin)*0.2,trmax+(trmax-trmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","tr","");
	cpgslw(7);
	cpgpt(nfile,xr,try,-1);
	cpgslw(2);
	cpgline(nfile,xr,try);
	//cpgerry(nfile,xr,try1,try2,1.0);
	cpgend();

	cpgbeg(0,"?",1,1);
	cpgpap(11.0,6.0/11.0);
	cpgslw(4);
	cpgsvp(0.1,0.95,0.7,0.95);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,logcn2min-(logcn2max-logcn2min)*0.2,logcn2max+(logcn2max-logcn2min)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","logC\\u2\\d\\dn","");
	cpgslw(7);
	cpgpt(nfile,xr,logcn2y,-1);
	cpgslw(2);
	cpgline(nfile,xr,logcn2y);
	//cpgerry(nfile,xr,logcn2y1,logcn2y2,1.0);

	cpgslw(4);
	cpgsvp(0.1,0.95,0.4,0.65);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,Vsmin-(Vsmax-Vsmin)*0.2,Vsmax+(Vsmax-Vsmin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	cpglab("","Velocity(km/s)","");
	cpgslw(7);
	cpgpt(nfile,xr,Vsy,-1);
	cpgslw(2);
	cpgline(nfile,xr,Vsy);
	//cpgerry(nfile,xr,Vsy1,Vsy2,1.0);

	cpgslw(4);
	cpgsvp(0.1,0.95,0.1,0.35);
	cpgswin(-xr[nfile-1]*0.05,xr[nfile-1]+xr[nfile-1]*0.05,Smin-(Smax-Smin)*0.2,Smax+(Smax-Smin)*0.2);
	cpgbox("bcnst",0.,0.,"bcnst",0.0,0.);
	char xaxis[64];
	sprintf(xaxis,"MJD-%.2f",MJD[0]);
	cpglab(xaxis,"Flux Density(mJy)","");
	cpgslw(7);
	cpgpt(nfile,xr,Sy,-1);
	cpgslw(2);
	cpgline(nfile,xr,Sy);
	cpgend();

}

float ReadS(char *filename,char *filename2)//Read flux den­sity function (filename(*.txt)，filename(*S.txt)）
{
	int i,j,nsub,nchan,nbin;
	float bw,length,tsub,dm,mjd,Dist;

	char file1[5],sub1[4],chan1[5],pol1[4],MJD1[4],freq1[5],smJy1[7],w101[4],w501[4],sn1[5],filename11[40],name[11];
	char file[64];
	int sub,chan,pol;
	float MJD,freq,smjy,w10,w50,sn,S=0.0;
	FILE *fin=fopen(filename,"r");
	if(fin==NULL)
	{
		printf("Can not open file %s,please check!\n",filename);
		exit(1);
	}
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd,&Dist);
	fclose(fin);

	FILE *fp1 = fopen(filename2,"r");
	if (fp1 == NULL)
	{ 
		printf("Can not open file %s,please check!\n",filename2);
		exit(1);

	}

	fscanf(fp1,"%s %s %s %s %s %s %s %s %s %s\n",file1,sub1,chan1,pol1,MJD1,freq1,smJy1,w101,w501,sn1);
	for (i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
			fscanf(fp1,"%s %d %d %d %f %f %f %f %f %f\n",file,&sub,&chan,&pol,&MJD,&freq,&smjy,&w10,&w50,&sn);
			S+=smjy;
		}
	}
	fclose(fp1);
	return(S);
}

float MIN(float x,float y)
{
	float z;
	z=x<y?x:y;
	return(z); 
}
float MAX(float x,float y)
{
	float z;
	z=x>y?x:y;
	return(z); 
}

int Gsacf(char *filename,float *c,float *SD1,float *v1,float *Dist1,float *MJD,char *name1)//Gaussian fits function（filename，coefficient(c1,c2,c3),fitting error,pulsar frequency,pulsar distance,MJD,pulsar name）
{
	int i,j,k,nsub,nchan,nbin;
	char filename11[64],name[11],title[64];
	float bw,freq,length,tsub,dm,mjd,Dist,MIN(float,float),MAX(float,float),ACF();

	FILE *fin=fopen(filename,"r");
	if(fin==NULL)
	{
		printf("Can not open file %s,please check!\n",filename);
		exit(1);
	}
	fscanf(fin,"%s %s %d %d %d %f %f %f %f %f %f %f\n",filename11,name,&nsub,&nchan,&nbin,&bw,&freq,&length,&tsub,&dm,&mjd,&Dist);
	float a[nsub],b[nchan],A[nsub][nchan],Agray[nsub][nchan],S[nsub][nchan],AMAX=0,AVERAGE=0,td0,vd0;
	int zeron=0;//the number of intensity value of 0
	for(i=0;i<nsub;i++)
	{
		for(j=0;j<nchan;j++)
		{
			fscanf(fin,"%f %f %f\n",&a[i],&b[j],&A[i][j]);
			AMAX=(AMAX>A[i][j]?AMAX:A[i][j]);//Find the intensity of the maximum
			AVERAGE+=A[i][j];
			if (A[i][j]==0)
			{
				zeron++;
			}
		}
	}
	AVERAGE=AVERAGE/(nsub*nchan-zeron);//The average intensity values
	//normalize intensity values  
	for (i=0;i<nsub;i++)
	{
		for (j=0;j<nchan;j++)
		{
			if(A[i][j]==0)
			{
				S[i][j]=0.0;
			}
			else S[i][j]=A[i][j]-AVERAGE;
		}
	}
	fclose(fin);
	
	float C[nsub][nchan*2-1],xr[nchan*2-1],kx[nchan*2-1],yr[nsub],ky[nsub],xrMIN=1000,yrMIN=1000;
	//normalized two-dimensional ACFs values,xr is x-axis range in the bottom of figure,kx is y-axis range in the bottom of figure,yr is x-axis range in the right of figure,ky is y-axis range in the right of figure
	for (i=0;i<nsub;i++)
	{
		for (j=-nchan+1;j<nchan;j++)
		{
			C[i][j+nchan-1]=ACF(j,i,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);//normalize ACF values
		}
	}
	for (j=0;j<nchan;j++)
	{
		kx[j]=(j-nchan+1)/(nchan/bw);
		kx[j+nchan-1]=j/(nchan/bw);
		xr[j+nchan-1]=ACF(j,0,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		xr[j]=ACF(nchan-j-1,0,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		xrMIN=MIN(xr[j],xrMIN);
	}
		
	for (j=0;j<nchan;j++)//Find half-width at half-maximum of the ACF along the frequency-lag axis at zero time lag
	{
		if (xr[j+nchan]<0.5)
		{
			vd0=j*(bw/nchan);
			break;
		}
	}
	for (i=0;i<nsub;i++)
	{
		ky[i]=i*(length/(nsub*60));
		yr[i]=ACF(0,i,nsub,nchan,S)/ACF(0,0,nsub,nchan,S);
		yrMIN=MIN(yr[i],yrMIN);
	}
	for (i=0;i<nsub;i++)//Find the time lag at zero frequency lag at the point where the ACF is 1/e of the maximum
	{
		if (yr[i]<1/e)
		{
			td0=i*(length/(nsub*60));
			break;
		}
	}

	int t=td0/2,v=vd0;//Set the fitting range
	if (t<4) t=4;//In order to avoid fitting too few data
	if (v<1) v=2;
	double Clog[t*2*v],D1[t*2*v][3],DD[3][t*2*v],U[t*2*v][t*2*v],V[3][3],X[3],Clog2[t*2*v];
	for (i=0,k=0;i<t;i++)
	{
		for (j=-v;j<v;j++)
		{
			if (C[i][j+nchan]<=0) Clog[k]=log(0.000001);
			else Clog[k]=log(C[i][j+nchan]);
			D1[k][0]=(bw/nchan*j)*(bw/nchan*j);
			D1[k][1]=(bw/nchan*j)*(i*(length/(nsub*60)));
			D1[k][2]=(i*length/(nsub*60))*(i*length/(nsub*60));
			k++;
		}
	}

	ginv(D1,t*2*v,3,DD,0.0001,U,V,t*2*v+1);//call function of calculate generalized inverse of a matrix 
	trmul(DD,Clog,3,t*2*v,1,X);//call matrix multiplication function
	trmul(D1,X,t*2*v,3,1,Clog2);

	double SD=0.0;//fitting error
	for (i=0;i<t*2*v;i++)
	{
		SD+=(Clog2[i]-Clog[i])*(Clog2[i]-Clog[i]);
	}
	*SD1=sqrt(SD/(t*2*v))/(t*2*v);
	*v1=freq;
	*MJD=mjd;
	*Dist1=Dist;
	for (i=0;i<strlen(name)+1;i++)
	{
		name1[i]=name[i];
	}
	for (i=0;i<3;i++)
	{
		c[i]=-X[i];
	}
}

float ACF(int dv,int dt,int nsub,int nchan,float *S)//calculate ACF value（Frequency shift，Time shift，number of subintegration，number of frequency channel，normalized matrix）
{
	float acf=0;
	int i,j,N,zeron=0;
	N=(nchan-abs(dv))*(nsub-dt);//the number of correlated pixel pairs
	if (dv>=0)
	{
		for (i=dt;i<nsub;i++)
		{
			for (j=dv;j<nchan;j++)
			{
				if(S[i*nchan+j]==0 || S[(i-dt)*nchan+j-dv]==0)
				{
					zeron++;
				}
				acf+=S[i*nchan+j]*S[(i-dt)*nchan+j-dv];
			}
		}
	}
	if (dv<0)
	{
		for (i=dt;i<nsub;i++)
		{
			for (j=0;j<(nchan-abs(dv));j++)
			{
				if(S[i*nchan+j]==0 || S[(i-dt)*nchan+j-dv]==0)
				{
					zeron++;
				}
				acf+=S[i*nchan+j]*S[(i-dt)*nchan+j-dv];
			}
		}
	}
	if (N-zeron>0)	acf=acf/(N-zeron);
	
	return(acf);
}
